The role of winding numbers in quantum Monte Carlo simulations 
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We discuss the effects of fixing the winding number in quantum Monte Carlo simulations. We 
present a simple geometrical argument as well as strong numerical evidence that one can obtain 
exact ground state results for periodic boundary conditions without changing the winding number. 
However, for very small systems the temperature has to be considerably lower than in simulations 
with fluctuating winding numbers. The relative deviation of a calculated observable from the exact 
ground state result typically scales as T 1 , where the exponent 7 is model and observable dependent 
and the prefactor decreases with increasing system size. Analytic results for a quantum rotor model 
further support our claim. 



One approach to numerical many-body physics is to 
stochastically sample the "world line" configurations of 
a real-space Euclidean path integral. The most cof. 
mon of these quantum Monte Carlo (QMC) methods 
are based on the Trotter, decomposition formula in dis- 
cretized imaginary time, til with a resulting systematic e&. 
ror that vanishes as the discretization At is decreased BeI 
For lattice models, methods have recently been caru 
structed for simulations in continuous imaginary time,ETEl 
hence directly giving results exact to within statistical 
fluctuations. A related method is the stochastic series 
expansion techniques (a generalization of Handscomb's 
methodty), which samples the power series expansion of 
the density matrix, and also is exact. For all these meth- 
ods, the configurations for periodic systems can be la- 
beled by a topological "winding number" w, which counts 
the net number of times the world lines wrap around the 
system in the course of propagation in imaginary time. 
In practice, it is often not possible to sample all winding 
number sectors, since changing w requires simultaneous 
modification of a number ~ L world lines, where L is 
the linear size of the system, with a resulting low accep- 
tance rate if L is large. Such noij^ecgodicity is clearly re- 
lated to the boundary conditiorJa-EJ (restricting to, e.g., 
w = can be considered a particular boundary condi- 
tion), and therefore results scaled to infinite system size 
are still correct, although there can be significant devia- 
tions from the exact periodic boundary condition results 
for any given small system size. The winding number is a 
consequence of the path integral formulation of quantum 
mechanics, and a fixed w is not related to the Hamilto- 
nian in a simple way. In many cases it is of interest to 
obtain approximation-free results for periodic boundary 
conditions, specifically, and the restriction to fixed w can 
then be a practical limitation of the QMC method. 

We here point out that the exact ground state can 
in fact be obtained in any winding number sector, even 
for the smallest possible systems. We base our claim on 
simple geometric considerations, and provide supporting 
simulation results for several many-body Hamiltonians 
(one- and two dimensional Heisenberg spin models) as 



well as analytic results for a quantum rotor. We expect 
our arguments to be generally valid for models for which 
the path integral is positive definite (i.e., there are no 
sign problemsQ'E3) , which is the case for, e.g., ID fermion 
systems and non-frustrated spin and boson systems in 
any dimension. 

In order to define the winding number, we first consider 
the standard path integral formulation of a continuous 
ID system. The partition function for a single particle 
of mass M in a potential V(x) ialj 



Z= I D[x(T)]e- s , 

>x(/3)=x(0) 



where the action is 



S= dr 
Jo 



M ( dx{r) 



2 V dr 



V(x(r)) 



(1) 



(2) 



The integral in (Q) is over all paths (or world lines) start- 
ing at position x at imaginary time r = and ending at 
the same position x at r = 0, where (3 denotes the inverse 
temperature. For a periodic system, the paths can be di- 
vided into topologically distinct classes, characterized by 
a winding number w defined as the net number of times 
the world line spatially wraps around the system. For 
a many-particle system, the winding number is defined 
as the net total displacement to the "left" or the "right" 
of the world lines in the course of propagation between 
t = and t = fl, divided by the length of the system. 
In higher dimensions, there is a winding number asso- 
ciated with each dimension, and the definition of these 
is a direct generalization of the above definition in one 
dimension. 

For interacting many-body systems, analogous real- 
space path-integrals (or related sumstl based on series ex- 
panding the density matrix operator e~@ H ) can be con- 
structed, and are suitable for numerical simulations in 
cases where the weight associated with the paths is pos- 
itive definite. Such QMC methods have been developed, 
e.g., for lattice fermions in one dimension!, latticeu-3 and 
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continuumH bosons, and quantum spin systems.^ In the 
case of fermions in higher dimensions, the path integral 
cannot be efficiently sampled directly, due to the non- 
positive definiteness of t he , w eight, which leads to the 
infamous "sign problem" £3li3 

Stochastic sampling of the bosonic paths within a sec- 
tor of a fixed winding number can be accomplished by 
local modifications of the world lines, and is typically 
a relatively straight-forward task. The global modifica- 
tions required in order to change the winding number are 
often practically impossible to carry out efficiently, how- 
ever. Recently "loop-cluster" algorithms have been de- 
veloped which in principle automatically sample all wind- 
ing number sectors.t3 However, such algorithms do not 
perform well for all models, and therefore the restriction 
to a sector with fixed winding number remains the only 
option in many cases. We note that the winding num- 
ber itself is related to long-range coherence in the system. 
For boson and spin systems, the winding number fluctua- 
tion (w ) is directly proportional toihe superfluid density 
and the spin stiffness, respectivelycl In some cases these 
quantities can, however, also be computed in a restricted 
simulation.ll-Ja 

We here argue that the exact ground state can actu- 
ally be studied in any w sector. Consider first again a ID 
periodic system with a single particle. A path integral 
configuration can be visualized as a string on the surface 
of a cylinder with periodic boundary conditions (i.e., a 
torus). As the temperature is lowered the length /3 = 1/T 
of the system in the imaginary time direction becomes 
much larger than the spatial system size, and the string 
can then wrap around the cylinder multiple times in both 
directions between r = and (3. In the winding number 
w sector, the net number of revolutions has to be w. Nev- 
ertheless, locally, in an imaginary-time segment of length 
At <C (3, it would not be possible to detect the effects of 
this restriction. Hence, any quantity that can be defined 
on a finite segment of the cylinder should be the same in 
any winding number sector as (3 — > oo. Since correlation 
functions decay exponentially with the imaginary-time 
separation in a finite system, due to the finite-size gap, 
all quantities which are not defined in terms of the global 
winding number itself should become exact as (3 — ► oo. 
This argument can clearly be generalized for a many- 
body system in any dimension, again of course provided 
that all paths have positive weights. In the case of mixed 
signs the positive and negative contributions to the par- 
tition function cancel as (3 — > oo, and physical quantities 
are given by finite ratios of two vanishing numbers. It 
is then not clear that all the winding number sectors be- 
come equal as (3 — * oo (although we have not proved the 
contrary), and we will not consider this intricate issue 
further here. 

We can rigorously prove that fixed w gives the correct 
ground state of a single particle in one dimension. This 
system (with ML 2 = 2ir 2 ) is equivalent to the quantum 
rotor, described by the Hamiltonian 



H = -d 2 /d0 2 . (3) 
The eigenfunctions labeled by the angular momentum m 



^ m (9) = exp(im9). 
The partition function is 

Z = E exp(-/W 



(4) 



(5) 



and can be transformed from a sum over angular mo- 
menta to a sum over winding numbers in the follow- 
ing way: In the discrete path integral formulation with 
At = (3/N the partition function is 



Z 



N 



DflH^^lexpf-Arff)!^) 



N 



do n 



j — l L Tfl 



cxp(— Arm 2 ) exp[im(0j+i — 9j)] 



(6) 



where At = (3/N and 0/v+i = 9\. Using Poisson's sum- 
mation formula, 



E /H= E F ( 27rn )> 



where F(n) is the Fourier transform of f(m), we obtain 

- 27m) 2 



N 

DO ft 



E ex P 



'3 + 1 
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Next we write the j (or t) dependence of 9 as 



6 j = mod[6 1 + —(j-l) + 



,2tt1 



(8) 



(9) 



where w is the winding number. For the time being we 
neglect the fluctuations 59 j, which are the same for all 
w. As At — * 0, only the term n = contributes in the 
above sum, unless 9 crosses the boundary at 9 = 0, in 
which case n = -lorn=+l will give the contributing 
term. Hence the partition function can be expressed as 



Z = .f(f3) E exp(-7rV/7?), 



(10) 



IV — — oo 



where f{(3) is a function containing the effects of the fluc- 
tuations 59 j which we have so far neglected. The easiest 
way to find f(/3) is to simply equate the above result 
with Eq. (|). This gives our final result for the partition 
function expressed as a sum over winding numbers: 



z = E V^7Pexp(^n 2 w 2 /(3). 



(11) 



W — — 00 
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This result could have been directly obtained using Pois- 
son's summation formula, Eq. ([?]), on Eq. (||). The above 
derivation shows that the new summation index indeed 
is the winding number, which hence can be thought of as 
a quantity conjugate to the angular momentum. 

The energy E = —{l/Z)dZ/dj3 given as a sum over 
angular momenta is 

E = -i^m 2 exp(-/?TO 2 ), (12) 

m 

and expressed as a sum over winding numbers 

w " " 

In Fig. H the energy of the rotor is plotted versus (3 on a 
log-log scale. It decreases as at high temperatures, 
changing to an exponential decay around (3=1. In the 
same graph we also show the energy evaluated in the w = 
sector. At high temperatures the results coincide with 
the full energy, but around /3 = 1 the behavior changes 
to be of the form /3 -3 / 2 , instead of exponential, as can 
be readily extracted from Eq. (|l3|). Hence, we indeed get 
the correct energy (namely, zero) in the w = sector, but 
the approach to the ground state is much slower than in 
the "ensemble" with fluctuating w. 
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FIG. 1. Rotor model : The full energy and the energy for 
the w = sector vs. j3. The lines are the assymptotic high 
and low (for w = 0) temperature forms. 

Next we present numerical results for several many- 
body models. We study the Heisenberg model, defined 
by the Hamiltonian 

# = ]TS,.S„ (14) 

where indicates a pair of nearest neighbors lattice 
sites, and S, is a spin-S* operator at site i. We consider 
the ID model with 5 = 1/2 and S = 1, as well as the 



2D model with S = 1/2. We note that these quantum 
spin models are formally equivalent to constrained boson 
systems. 

For the simulations we use the stochastic series ex- 
pansion algorithm^ which is based on a power series 
expansion of e~@ H , and hence is not a standard path 
integral method. The configuration space is neverthe^. 
less very strongly related to an Euclidean path integral,u 
and the winding number has exactly the same meaning. 
For systems of linear size L < 16 Monte Carlo updates 
changing the winding number can be easily carried out, 
but for larger systems a restriction to fixed w is necessary 
in practice. The advantage of the method is that there 
are no other sources of systematic errors. We present 
energy results both for w = and fluctuating w, and 
compare with exact diagonalization data. 

The simulation scheme is formulated in the standard 
basis where the operators Sf are diagonal. The internal 
energy per spin can be calculated in two different ways 
in the simulation; from the nearest-neighbor correlation 
function (S*S* +1 ) as well as from a manifestly rotation- 
ally invariant estimator giving the full (Si • Sj+i).u We 
define 

E, = 3(S?Sf +1 ), (15a) 
E s = ^ ■ S< + i). (15b) 

The agreement between the two estimates can serve as 
a good internal check of the spherical spin symmetry in 
the simulation. With w fixed in a finite system, the spin- 
rotational symmetry is broken since only the xy-terms 
are involved in the spatial propagation (spin flipping) 
of the path and the estimates will therefore not agree 
completely. 

In Fig. E s and E z are graphed versus the inverse 
temperature (3 used in simulations of an 8-site S = 1/2 
chain, both in the w = sector and with fluctuating 
winding numbers. In the w = sector, the rotation- 
ally invariant and diagonal estimates approach the exact 
ground state result from above, and below, respectively. 
For fluctuating winding numbers the exact ground state 
energy is obtained within statistical errors already for 
j3 » 8 for this system size, reflecting the large finite-size 
gap and the exponential approach to the ground state 
with increasing 0. 

In Fig. |3| the difference between Monte Carlo data 
obtained in the w = sector and the exact temperature- 
dependent energy is graphed on a log-log scale for L = 
4, 8 and 16. At high temperatures the difference tends to 
vanish rapidly since the system is too short in the imag- 
inary time dimension to allow for the winding number 
to change from zero. At low temperatures the devia- 
tion decreases as /3 -7 , with 7=1 within the numerical 
precision, and there is a maximum at an intermediate 
value (3* . As expected, the maximum difference decreases 
rapidly with system size. The relative error decreases 
roughly as l/L at a given temperature, and (3* increases 
roughly as L. (3* hence reflects the inverse finite-size gap. 
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FIG. 2. ID spin- 1/2 model: The energy estimators E 3 and 
E z vs. 13 for L=8. 



FIG. 4. ID spin-1 model: The deviation from the ground 
state energy of E s vs. /3 for L=4 and 8. The lines have slope 
-1. 



These results confirm that for a given desired accuracy 
of the energy at a given temperature, there will be some 
system size beyond which it is not necessary to sample 
different winding number sectors. Furthermore, the ex- 
act ground state is always obtained as [3 — > oo. However, 
for small systems we have to use higher values of (3 in 
order to achieve a certain accuracy if we do not sample 
all winding number sectors. 
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FIG. 3. ID spin-1/2 model: The deviation from the exact 
energy of E s vs. j3 for L=4,8 and 16. The lines have slope -1. 

Next we consider the ID S = 1 model and the 2D 
S = 1/2 model. The relative deviation of the w = 
energies from the exact ground state values are graphed 
versus [3 for different system sizes in Fig. | and Fig. 
H (for the 2D systems with L > 4 "exact" results were 
obtained in simulations with fluctuating w) . Again we see 
that at low temperatures the ground state is approached 
as instead of the exponential approach expected for 
the exact energy. 

To summarize we have given an intuitive argument why 
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FIG. 5. 2D spin-1/2 model: The deviation from the ground 
state energy of E s vs. /3 for L=4,8 and 16. The lines have 
slope -1 

restricting the winding number in QMC simulations num- 
ber should not affect calculated ground state properties 
of even the smallest lattices. We have given strong nu- 
merical evidence to support this statement, in addition to 
rigorous results for a simple quantum rotor model. Typ- 
ically, the asymptotic deviations from the exact periodic 
boundary condition results scale as /3 -7 at low temper- 
atures, with a prefactor that goes to zero as the system 
size increases. Hence, in terms of the (3 needed to obtain 
the ground state, to the statistical precision possible in 
QMC simulations, there does not appear to be any ad- 
vantage of fluctuating numbers for moderate and large 
systems. 

Our considerations are closely related to the 
imaginary-time boundafjf conditions recently discussed 
by Tauber and NelsonO They showed that the condi- 
tion that the world line configurations at r = and r = (3 
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being equal can be relaxed, still giving the same result 70, 875 (1993). 

as the periodic imaginary time boundary condition when 17 U. C. T auber and D. R. N elson, to appear in Phys. Rep. 
[3 — > oo. In analogy with our discussion above, the rea- (1997) ( |cond-mat/9608057| ) . 

son for this is that the changed boundary condition does 
not affect the behavior of the world lines in a segment of 
length <C /?. 
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